%% Codes for the article "Beliefs and Financial Crises", by Arvind Krishnamurthy and Wenhao Li.
% This file generates the parameter combinations using the Smolayak grid.
% Last edit: Feb 2, 2020
function   [index_minimum_error,parameters_matrix,dis_vec,round_index_best] = calibration_Smolayak_grid_method(  model_option, starting_round, ending_round, output_best_version )
% if output_best_version=true, then simply output the best version so far from the simulation results. 

if(nargin<2)
    starting_round = 1;    
end
if(nargin<3)
    ending_round = starting_round;    
end
if(nargin<4)
    output_best_version = false;
end

%% Choose whether to use the behavioral or rational model
model = model_option.model_name;
disp( [ '-----------------------------------------------------------------------' ]  )
disp( [ '                                 Calibration of the ',  model  ,' model  ' ]  )
disp( [ '-----------------------------------------------------------------------' ]  )


%%  Initialization: generate default parameters, and a row of parameters for estimation
target_moment_file = 'moment_target.xlsx'; 
target_values =  readtable( target_moment_file, 'Sheet',1,  'Range',  ['B3:B9'] ,'ReadVariableNames',false).Variables;
bank_leverage_target = target_values(7);
deviation_percentage_vec = [ 0.1, 0.05, 0.025 ];    % important: the range of parameter values

if( model=="benchmark" )
    shutoff_positions = [4,5,7];   % for the benchmark model, only the average frequency of the liquidity shock matters. Thus, we shut off position 4 for lambda L to H and Position 5 for lambda H to L.
    % Note: we shut off position 7, eta, which is adjusted afterwards to match average leverage. 
    accuracy_level=2;      
elseif( model=="rational" )
    shutoff_positions = [7];  
    accuracy_level=2;      
else
    shutoff_positions = [7];    
    accuracy_level=2;      
end

if(isfield(  model_option, 'shutoff_positions'  ))    % note: if shut off
    shutoff_positions = model_option.shutoff_positions;   
end

    
%% ********** Start the calibration for the model.  In total we have rounds that are specified by the length of "deviation_percentage_vec". 
for( round_index = starting_round:ending_round )  % shrink the parameter deviation from the optimal parameter combination form the last round.

if(output_best_version && round_index>starting_round)  % if output the best version, then do not execute the next round.  Only excute the current round
    break;
end

%% Step 0: Generate parameter matrix
% We find the baseline parameter for the iteration in two ways: (1) use the current best results.  (2) Use a default
% results.  We use (1) whenever possible. 
if( round_index==1 ) 
    % note: the first round should start from initialization
    S = load( ['solutions/initialization/',model ,'_init.mat'] );
    par_baseline = S.model_sol.par;
else
    % future rounds should use the best result from previous rounds.
    S = load( ['solutions/baseline_',  model ,'.mat'] );
    par_baseline = S.model_sol.par;
end
par_baseline.alpha = 0.05;    % alpha is a parameter directly set by the calibration.
par_baseline.model_name = model;   
par_baseline = initialization(par_baseline);    % Need initialize again with the added model name and crisis frequency.
par_baseline.alg_control.solve_everything_again = isfield( model_option, 'solve_credit_spread_again' ) && model_option.solve_credit_spread_again ;
parameters_baseline = [ par_baseline.AL, par_baseline.AH - par_baseline.AL ,  par_baseline.lambda_H,  par_baseline.lambdaLH,	par_baseline.lambdaHL,  par_baseline.sigmaK,	par_baseline.eta];
% In total there are 7 parameters to be calibrated. 
% Note: parameters_row needs to be consistent with the "update_par_with_vector" function. 
main_file_name= strcat('calibration files/calibration_smolayak_', model ,'.xlsx');
parameters_row = parameters_baseline;   
deviation_percentage= deviation_percentage_vec(round_index); 
parameters_matrix = smolayak_parameter_grid_generation( parameters_row,  accuracy_level, deviation_percentage,  shutoff_positions );
N_row = size(parameters_matrix,1);

% Important: If the corresponding excel file has already parameters in it, we should not change it. Instead we just use
% whatever is already in the file. This avoids mistakes. 
start_row = 5;   end_row = start_row+N_row-1;
parameter_input = readtable( main_file_name , 'Sheet', round_index,  'Range',   [ 'D', num2str(start_row), ':J', num2str(end_row) ], 'ReadVariableNames', false  );
parameter_input = table2array(parameter_input);
if( sum(sum(~isnan( parameter_input ))) > 0 )   % As long as there are data in the table, simply use the original table.
    parameters_matrix = parameter_input;
end

% Then write the matrix to excel.  If we already have content in the excel, then this operation does not change the
% content.
parameters_matrix_output = [   (1:N_row)' , parameters_matrix  ];  
writetable( array2table(parameters_matrix_output) , main_file_name , 'Sheet',round_index,'Range',   [ 'C', num2str(start_row), ':K', num2str(end_row) ]  ,'WriteVariableNames',false)

non_na_rows = sum( isnan(parameters_matrix)' )'==0;
parameters_matrix = parameters_matrix(non_na_rows,:);


%% Step 1: Solve and Simulate Different Scenarios
if(~output_best_version)

solution_path = strcat('solutions/smolayak_', model ,'/round', num2str(round_index), '/') ; 

% (2) Solve the model for each parameter combination.
p_init = S.model_sol.p_matrix;    
solve_model_again =  isfield( model_option, 'solve_model_again' ) && model_option.solve_model_again; % If we want to solve again regardless
smolayak_solve_parameter_combinations( parameters_matrix,  par_baseline,  p_init, solution_path, solve_model_again)

% (3) Adjust eta and kappa0 so that the stationary distribution of the model implies  the bank levearge ratio the same as
% moment target. 
iter_end = size(parameters_matrix,1);  %N_trials
iter_start=1;
if( (~isfield( model_option, 'not_adjust_eta' )) || (~model_option.not_adjust_eta) )  % This section can be jumped if model_option has the not_adjust_eta field or it is false. 
    parsave = @(file_name, model_sol) save(file_name, 'model_sol');
    loss_given_default = par_baseline.loss_given_default;
    warning('off','all');
    for(iter_combination=iter_start:iter_end)   
        file_name = [solution_path, num2str(iter_combination) ,'.mat'];
        disp( '-------------------------------------------------------' )
        disp( ['       Adjust eta and kappa0 for parameter combination   ', num2str(iter_combination)] )
        disp( '-------------------------------------------------------' )
        S = load(file_name);   model_sol = S.model_sol;      model_sol.par.model_name = model_option.model_name;
        try
            [eta, hat_kappa0, lvg, kappap_multiplier] = adjust_eta_and_kappa0( model_sol, bank_leverage_target, loss_given_default );
            disp( ['Final adjusted leverage=', num2str(lvg)] ); 
            model_sol.par.eta = eta;
            model_sol.par.hat_kappa0 = hat_kappa0;
            model_sol.par.kappap_multiplier = kappap_multiplier; 
        end
        disp( ['**** adjusted eta=', num2str(eta,2), ', hat_kappa0=', num2str(hat_kappa0,2), ',  for "', file_name, '"'] );
        parsave(file_name, model_sol);
    end
end

% Then output the eta values to the excel  .  
eta_vec = ones(iter_end,1);    iter_start=1;
for( iter_combination=iter_start:iter_end )
    file_name = [solution_path, num2str(iter_combination) ,'.mat'];
    S = load(file_name);
    eta_vec(iter_combination) = S.model_sol.par.eta;
end
writematrix( eta_vec , main_file_name,  'Sheet',  round_index ,  'Range',  ['J', num2str(4+1) ,':', 'J', num2str(4+iter_end) ]  );  

% (4) Solve for credit spread
parfor( iter_combination = 1:iter_end )
    file_path = [solution_path, num2str(iter_combination), '.mat'];
    S = load(file_path);       par = S.model_sol.par;   par.alg_control.solve_everything_again = par_baseline.alg_control.solve_everything_again;
    par.alg_control.show_progress=true;  par.alg_control.plot=false;
    if(isfile(file_path))  % if the file exists, try to add credit spread onto it.
        try  
            credit_spread = post_processing_define_spreads( file_path, par );   % After saving the baseline version, then add credit spread to it, and overwrite the file
            solved = sum(sum(credit_spread ~= real(credit_spread)))==0;   % Whether credit spread solved has all elements being real rather than complex numbers. If the latter, then the algorithm didnt' succeed.
            if(~solved)
                disp( [ file_path, ' cannot be solved for credit spread!' ] )
            end
        catch
            disp( [ file_path, ' cannot be solved for credit spread!' ] )
            solved = false;
        end
        if( ~solved )
            try
                par.use_smoothing_for_v = true;
                post_processing_define_spreads( file_path, par ); 
            catch
                disp( [ file_path, ' cannot be solved for credit spread, using the smoothing method!' ] )
            end
        end
    end
end


if(par_baseline.alg_control.solve_everything_again)
for(  iter_combination = 1:iter_end )
    file_path = [solution_path, num2str(iter_combination), '.mat'];
    try
        S = load(file_path); 
        path_cs_plot = [ strcat('solutions/smolayak_', model ,'/round', num2str(round_index), '_credit_spread/') ];
        mkdir(  path_cs_plot );
        h = figure( 'Position', [0 0 600 400] ) ; 
        subplot(1,2,1);plot( S.model_sol.par.w_grid, S.model_sol.credit_spread );  xlabel('w'); ylabel('credit spread');
        subplot(1,2,2);plot( S.model_sol.par.lambda_grid, S.model_sol.credit_spread );   xlabel('\lambda'); ylabel('credit spread');
        saveTightFigure( h,  [path_cs_plot, num2str(iter_combination), '.pdf' ] );
        close all;
    catch
        disp( [file_path ' credit spread cannot be plotted. ' ])
    end
end
end

% (5) Simulate different combinations of parameters. 
output_info.file_name =  main_file_name;  output_info.sheet=round_index;  
output_info.starting_cell_column = 'M';     output_info.starting_cell_row = 5;  
output_info.ending_cell_column = 'AF';
smolayak_simulate_parameter_combinations( parameters_matrix, solution_path, output_info );

end  % end of Step 1.


%% Step 2: Find the parameter combination that generates the minimum calibration error 
n_row = size(parameters_matrix,1);
if( model=="benchmark" )
    column_output = 'C';
elseif( model == "rational")
    column_output = 'D';
else
    column_output = 'E';
end
moment_index = [ 1:7 ];
weights = ones(size(moment_index)); 
weights(1) = 3; % extra weight on the average liquidity premium. 
[index_minimum_error, dis_vec, round_index_best, calibrated_par_values] = smolayak_find_best_parameter_combination( model_option, 1:round_index,  n_row, weights );
display( ['Best solution in ', model, ' is round ', num2str(round_index_best), ' solution ', num2str(index_minimum_error), '.mat' ]  );

%% Step 3: Update the new baseline parameters and output the best version of the model
if(round_index_best<Inf)
    file_path = ['solutions/smolayak_', model ,'/round', num2str(round_index_best), '/',  num2str(index_minimum_error), '.mat'];
    S = load(file_path); 
    simulation_key_moments_for_parameter_calibration;
    m = simulation_results_output(S, 15000);  % the second parameter is simulation year
    model_sol = S.model_sol;
    % Save the updated baseline model
    save( [ 'solutions/baseline_', model, '.mat' ],    'model_sol'  );
     % Then write calibrated moments to the calibration moment file. 
    calibrated_moment_values_full = calibration_moments(m,S.model_sol.par);
    writetable( array2table(calibrated_moment_values_full([1:8,15:16])') , target_moment_file , 'Sheet', 1,  'Range',   [ column_output, '3:',  column_output, num2str(3+length(target_values)+1) ]  ,'WriteVariableNames',false);
    % Write the parameter values to the target_moment_file
    par = model_sol.par;
    par_values = [  par.AL, par.AH-par.AL,  par.lambda_H, par.lambdaLH, par.lambdaHL, par.sigmaK, par.eta   ];
    writetable( array2table(par_values') , target_moment_file , 'Sheet', 2,  'Range',   [ column_output, '16:',  column_output, num2str(16+length(par_values)-1) ]  ,'WriteVariableNames',false);
end

end